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FINITE ELEMENT MODELLING 
OF THREE-DIMENSIONAL STRUCTURES 


By 

William R. Case, Jr. 
Goddard Space Flight Center 


SUMMARY 

A three-dimensional finite element originally devised by A. Hrennikoff has 
been used to determine its accuracy in predicting the deflections, stresses and 
strains of an elastic continuum. This finite element model is of rectangular 
parallelepiped shape and is essentially a constant stress-strain element. In 
order to gauge its accuracy, a classic problem in three-dimensional elasticity 
with a known closed form solution was selected and this structure was modelled 
with the Hrennikoff finite element. The problem selected for this initial investi- 
gation was that of a prismatical bar stretching due to the action of its own 
weight. The results of the finite element analysis show excellent agreement 
with the theory even for a very coarse mesh size subdivision of the prismatical 
bar thereby indicating that the Hrennikoff finite element model can be used 
with confidence to model three -dimensionally elastic structures. Limitations 
on the application of this element are indicated in the text. 
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FINITE ELEMENT MODELLING 
OF THREE-DIMENSIONAL STRUCTURES 

William R. Case, Jr. 

Goddard Space Flight Center 


INTRODUCTION 

Finite element modelling techniques for the determination of displacements 
and stresses in elastic bodies have received considerable attention in the litera- 
ture throughout the last several years due to the development of acceptable 
numerical methods for computer solution of large order matrix equations. In 
particular, a great deal of effort has been concerned with the development of 
the stiffness characteristics of what may be termed one and two-dimensional 
elastic elements such as beams and plates while little has been done with three- 
dimensional elasticity problems. In recent years, however, the analyst has 
developed a need for obtaining solutions to complicated structures such as heat 
shields, rocket engines, solid propellant boosters, etc., which are truly three- 
dimensionally elastic in their load carrying capability. 

It is the intent of this project to investigate ways of modelling three- 
dimensional structures in an effort to obtain a suitable finite element stiffness 
matrix which may be used in the NASA NASTRAN program for finite element 
analysis of structures. This report is concerned with the first method investi- 
gated and with the comparison of the results with theoretical solutions for one 
particular application. 

The three-dimensional finite element model selected for this first investi- 
gation was originally developed by A. Hrennikoff (Reference 1) and was selected 
for several reasons. It is one of the simplest approaches to obtaining a stiffness 
matrix to represent a three-dimensional elastic continuum, and also the cor- 
responding approach for two-dimensional elements (e.g. plate bending elements) 
has been shown to yield excellent results. 


FINITE ELEMENT MODEL 

The method developed by Hrennikoff for devising a finite element model of 
a three-dimensional elastic continuum consists of replacing this continuum by a 
framework of rods capable of carrying tension and compression and pin-joined 
at the corners as shown in Figure 1. By a judicious choice for the areas of these 
tension-compression members the framework can be made to deform according 
to Hooke’s law when it is loaded at its corners with concentrated forces that are 
statically equivalent to the uniform stresses on the elastic solid. 
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Figure 1. 3-D Hrennikoff Model 


The actual calculation of the areas of the rods of the framework will not be 
carried out here since it is covered in detail for two-dimensional frameworks 
in References 1 and 2. In general, however, the procedure is as follows. We 
load the framework with concentrated loads at the grid points which are stati- 
cally equivalent to the following uniform stress states: 

1. Uniform normal stress t in the x i direction and vo in the x j and 
x k directions. 

2. Uniform normal stress o in the x j direction and va in the x k and 
x A directions. 

3. Uniform normal stress a in the x k direction and v a in the x | and 
Xj directions. 

4. Uniform shear stress / in each of the planes x t - x. , Xj -x k and 

x k“ x i- 

From these will result nine independent equations for elongations or shear of the 
framework in the various directions which are equated to those prescribed by 
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Hooke's law for an elastic volume under the action of the same stresses. In- 
volved in these equations are nine parameters, namely the areas of the rods on 
the edges of the framework (A i} A. and A k ) parallel to the three axes; the areas 
of the rods that form the diagonals on the faces of the framework (A. , A. k and 
\ i ); the ratios of two of the lengths of the framework to the other length and 
Poisson's ratio. It is found from these nine equations that the length ratios are 
arbitrary (i.e. we can make the framework the same size as the continuum we 
are modelling) and that Poisson's ratio (of the elastic continuum we are attempt- 
ing to model) must have a value of: v 1/4. The resulting rod areas are found 
to be: 
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and a 0 is any reference length. The rather surprising result for the prescribed 
value of 1/4 for Poisson's ratio was originally discovered by Hrennikoff and is 
due to the fact that there are more stress conditions to satisfy than there are 
unknown rod areas. This difficulty can be avoided by adding more rods to the 
framework or by making the edge members beams that are rigidly connected 
to one another. For the purposes of this analysis, however, the value of =1/4 
is acceptable and further complication of the model in order to free the restric- 
tion on Poisson's ratio is not warranted at this time. 

If the rods of the framework have areas as given above, then the framework 
will deform exactly the same as the elastic solid if the solid is acted upon by 
uniform stresses and the framework is loaded at the grid points with statically 
equivalent concentrated forces. A three-dimensional elastic continuum can 
therefore be subdivided into small rectangular prisms and each of these in turn 
replaced by an equivalent framework whose rods have the areas prescribed above 
with the a i , a j( and a k the size of the frameworks (i.e. the distances between 
grid points). This equivalent structure, consisting of all of the frameworks, is 
loaded at its grid points with concentrated forces that are obtained by lumping 
the actual body and surface forces acting on the elastic continuum. If the sub- 
division of the body is fine enough, the deflections of the finite element model 
will approach the true deflections of the elastic body. 

Once the bar areas are known, the stiffness matrix of a single framework 
can be calculated by the standard methods. The development of the stiffness 
matrix is carried out in Appendix A for the 24 degrees of freedom (three trans- 
lation degrees of freedom at each of the eight grid points of the framework). 

This stiffness matrix could be programmed into the NASTRAN program so that 
the user would have only to call for this element and from the grid point geometry 
and material property data input by the user the program would automatically 
calculate the stiffness matrices for all of the elements (frameworks) that com- 
prise the structure. Since the general stiffness matrix, as shown in Table 1, 
has all of its coefficients in the form of simple algebraic equations, the calcula- 
tion of this matrix by the program would require only a relatively short time. 

Other three-dimensional finite elements may be more elaborate and for a 
given mesh size subdivision of a three-dimensional structure may give better 
results but would require more computer time to formulate the individual ele- 
ment stiffness matrices. Since the total computer time to solve a structures 
problem depends on the time to formulate the stiffness matrix for each element 
as well as the total number of grid points, it may be advantageous to use simple 
finite elements with a fine mesh subdivision of the structure as opposed to a 
coarser mesh with more elegant finite elements. At any rate, the simplicity of 
the Hrennikoff element, as depicted by its stiffness matrix in Table 1, is con- 
sidered a definite advantage at this point. 
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As mentioned above, the Hrennikoff element stiffness matrix could be auto- 
mated in the NASTRAN program although this has not been done and will not be 
done until other methods of modelling three-dimensional structures have been 
tried. In order to use the Hrennikoff element in the NASTRAN program at the 
present time the user must calculate the areas of all of the individual rods and 
input these to the program. The program does have rod elements and it will 
assemble all of the rods that make up the structure and solve the problem. 

TEST PROBLEM 

In general there are very few closed form solutions to three-dimensional 
elasticity problems that could be used as checkout problems for the three- 
dimensional finite element. Solutions that do exist fall into three major cate- 
gories: closed form solutions to three-dimensional problems with finite bound- 
aries, with boundaries at infinity, and approximate solutions to some classes of 
problems in the form of infinite series. Those in the first category are simple 
elasticity problems such as the bending of a three-dimensional beam, torsion of 
circular shafts and stretching of prismatical bars under the action of gravity 
loads. Within the second category lie such problems as the Boussinesq, Kelvin 
and Cerruti problems of forces acting on infinite or semi-infinite spaces. These 
latter problems can be used to check the finite element by making the model have 
specific finite dimensions and by applying conditions at the boundaries of the 
model that are the same as the exact elasticity solution prescribes at the cor- 
responding points of the infinite medium. Such a method is discussed in Refer- 
ence 3 in general and also specifically for the bending of a thick plate. 

As a first trial for the 3-D Hrennikoff model, it was decided to use one of 
the simplier elasticity problems, in particular the stretching of a prismatical 
bar under the action of gravity loads. This problem is completely described in 
Reference 4 and only the solution will be presented here. For a bar as shown 
in Figure 2, under the action of its own weight and supported at the top such that 
the vertical deflection is zero at the point x = y-0, z=C and with a stress 
distribution of: 



Figure 2. Prismatical Bar 


The deflections of the bar are: 

v y 

u = - T xz 

v y 

V = - T yz 

W = - [(z 2 - P 2 ) + v(x 2 + y 2 )] 

where y is the weight per unit volume of the bar. 

For convenience we will nondimensionalize the deflections and stresses to 

fyn 2 ) 

the maximum vertical deflection along the z axis \~2E~7 anc * the maximum 
normal stress in the z direction (y£): 
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Where 
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Equations (5) - (8) for the deflections and stresses are the exact solution to 
the bar hanging by its own weight and will be used as the basis of comparison for 
the finite element model results. The finite element model of the bar in Figure 2 
is obtained by selecting a mesh size and replacing the elastic continuum within 
each mesh by the Hrennikoff frameworks as shown in Figure 3. Note that since 
the theoretical solution for w is symmetrical with respect to x and y and u 
and v are antisymmetrical, we can model only one quarter of the bar as indicated 
in Figure 3. 

For the 1/4 model shown, a mesh size of 3 by 3 by 5 is indicated which means 
that the original continuum is now represented by 45 smaller volumes (3x3x5) 
each of which is a Hrennikoff framework of the type in Figure 1. The mesh size 
shown above was arbitrarily chosen and was used for several checkout problems 
in which the length of the bar was varied to determine whether the results would 
be more or less accurate for short stubby bars (in which Poisson's effect is large) 
than long slender ones. Other mesh sizes were chosen for a particular length bar 
to determine the rate of convergence of the finite element model results with mesh 
refinement. All computer runs were made on the IBM 7094/7040 DCS using the 
NASA NASTRAN program for finite element structural analysis by the displace- 
ment method (the force method will also be available at a later date). Solution 
times were all on the order of 10 minutes. The results of the finite element 
analyses together with the theoretical deflections and stresses for three different 
aspect ratio (a/P ) bars are shown in Figures 4-a through 6-c. All bars had a 
semi -width of a = 3", and the three aspect ratios were obtained from bars of 
length 50", 1.0", and 6". In all three cases, the finite element mesh size was 3 by 
3 by 5 as shown in Figure 3. The gravity loading, 

Z - y 


7 


z, w 



Figure 3. 


Finite Element Model of One Quarter of the Prismatical Bar (Mesh Size 3*3 X 5) 


acting on the continuum was simulated on the finite element model by integrating 
over a volume around each grid point in the mesh and lumping the resultant load 
at the grid point. For each volume in the mesh, the equivalent Hrennikoff frame- 
work was made up from pin-jointed rods with areas given by Equation (2), taking 
i, j , and k to be the x, y, and z directions, respectively. Note that these 
equations give the areas of exterior rods, that is, those between two grid points 
not shared by another volume in the mesh. For two grid points that share rods 
from more than one framework, only one rod is used between these grid points 
that has the total area of all rods from the frameworks that connect between 
these grid points. Thus, along the x direction, for example, two grid points 
that are interior to the solid we are modelling will have four frameworks that 
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Figure 4a. Vertical Deflections fora/2 - .06 (Finite Element Mesh Size 3x3x5) 
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Nondimensional Vertical Deflection 


Figure 5a. Vertical Deflections for a/£ — .3 (Finite Element Mesh Size 3*3*5) 
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Nondimensional Lateral Deflections 



Figure 5b. Lateral Deflections for a/£ - .3 (Finite Element Mesh Size 3x3x5) 
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Nondimensional Vertical Deflection 


Theoretical (x y 0) 



Figure 6a. Vertical Deflections for a/£ - .5 (Finite Element Mesh Size 3x3x5) 
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come together along a common edge thereby requiring that A. , as given in 
Equation (2), be increased by a factor of four to obtain the single rod area be- 
tween these two grid points. 

As indicated in Figures 4-a, 5-a, and G-a, the finite element results are in 
excellent agreement with the theory for the vertical deflections. Shown in these 
figures are the w deflections along the z axis as well as those along the furthest 
edge; x = y = 1 (except for the a /t - .06 case where the difference in the two 
is too small to show up). The w deflections along the x - y = 1 edge are where 
the maximum errors in any vertical deflection occur. Figures 4-b, 5-b, and 6-b 
show the nondimensionalized lateral deflection in the x direction and again the 
agreement is excellent except on the boundaries at z = 0 and 1.0 where the 
error increases linearly with x. Curves for the nondimensionalized lateral 
deflection in the y direction are not shown although they are exactly the same 
as those for the x direction if x is replaced by y in Figures 4-b, 5-b, and 6-b. 
That is, the finite element results yielded perfectly the symmetry condition: 

u(x, y, z) = v(y, x, z) 

which is due to the fact that 3u/3x and dv/dy arc identical as indicated by 
Equations (5) and (6). 

The results for the stresses in the bar are shown on Figures 4-c, 5-c, and 
6-c but were not obtained directly from the finite element computer runs. These 
were obtained by first calculating average strains using the grid point deflections 
and then converting these to stresses using Hooke's law. As such, they represent 
average stresses of the continuum over the volume bounded by the eight corner 
grid points of each framework. The figures show that the nondimensionalized 
axial stress, a z is in excellent agreement with the theory. All of the other 
stresses should be zero as indicated by Equation (4) and the results obtained 
for these from the grid point deflections indicated that the maximum stress 
other than c? z was always less than 2% of the maximum a z which is unity. 

For the bar with aspect ratio a/£ = .5, several finite element analyses were 
performed for various mesh sizes to investigate the rate of convergence of the 
results with model mesh refinement. Besides the 3x3x5 mesh discussed above, 
four other mesh sizes were employed, namely 1x1x2, 2x2x4, 3x3x6, and 
4x4x8. Therefore, if M is designated as the mesh size (the total number of 
bays that the bar is subdivided into) and N is designated as the number of sub- 
divisions along the x or y axis, then for these four problems: 

M = (N) (N) (2N) = 2N 3 
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where N =1,2, 3, and 4. For the case N = 4 (i.e., the 4x4x8 mesh) the re- 
sults for the vertical and lateral deflections and the axial Stresses are shown in 
Figures 7-a, 7-b, and 7-c. As before for the 3x3x5 mesh run of the bar with 
a/f = .5 , the finite element results are in excellent agreement with the theory 
with the boundary errors decreasing over those for the 3x3x5 mesh as would 
be expected. The solution time on the computer for the 4x4x8 mesh model was 
approximately 16 minutes as compared to 8 minutes for the 3x3x5 mesh. 

Figure 8 shows the comparison of several of the finite element results for 
the a /(. = .5 bar for the four mesh sizes selected (N = 1, 2, 3, and 4 in Equa- 
tion (10)). The vertical deflection at the origin, the results of which are desig- 
nated by the circled data points on Figure 8, are exact for all mesh sizes. This 
is due to the assumptions that went into the Hrennikoff model and the linearly 
varying stress state in the three-dimensional solid that is being modelled. 

The vertical deflection at the corner x = y = 1 and z = 0 is also shown 
since the largest magnitude error in any vertical deflection occurs here. As 
indicated by the curve, the finite element result converges to within approximately 
1% of the theoretical value for the finest mesh size analyzed. Similarly, the 
boundary deflection in the lateral direction at x = z = 1 (for all y ) converges 
rapidly to the theoretical value also. The strain energy curve shown on the 
figure was obtained from the grid point deflections and the input grid point forces 
and is included since it gives a sort of weighted or average error of the deflec- 
tions when compared to the theoretical strain energy. Even for the very coarse 
mesh size of 1 xlx2, all of the results shown on Figure 8 are much better than 
one would have anticipated. This is probably due to the nature of the problem 
investigated. Since the load is a distributed one which varies uniformly over 
the bar, one would expect to get better results for the same mesh size for this 
type problem than for one in which a single concentrated load or a nonuniform 
distributed load was acting. 


CONCLUSIONS 

On the basis of this single investigation, it is concluded that the Hrennikoff 
model may be used as a practical model for three-dimensional elasticity prob- 
lems in which the stress distributions and boundary conditions do not vary 
rapidly over the boay. It has been shown that for a bar hanging by its own weight, 
this technique for modelling three-dimensional structures gives excellent re- 
sults even for very coarse mesh sizes. Since the Hrennikoff model is a constant 
stress strain finite element, fine r mesh sizes will have to be employed for three- 
dimensional elasticity problems in which the stress gradients are larger than 
those investigated here. This disadvantage may well be offset by the distinct 
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Nondimensional Vertical Deflection 



Figure 7a. Vertical Deflections for a/l = .5 (Finite Element Mesh Size 4x4x8) 
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Nondimensional Lateral Deflections 






Ratio of Actual to Theoretical Values 



Note: Data points are for models with N bays in the x and y directions and 2N bays in the 
z direction (N 1, 2, 3 and 4) so that M 2N^. 


Figure 8. Comparison of Results for Various Mesh Sizes for Bar with a/e = 0.5 
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advantage of the simplicity of the Hrennikoff model. As indicated on Table 1, 
the stiffness matrix for any aspect ratio Hrennikoff model is available in a 
closed form solution involving only algebraic expressions. As more complexity 
is built into other three-dimensional finite elements, their corresponding stiff- 
ness matrices will be more complex and thereby require more computer time 
to formulate. Although the more complex finite element will allow problems to 
be solved with coarser mesh sizes, this savings in computer time may be offset 
by the increase in computer time required to formulate the stiffness matrix. 

Future work in this area of three-dimensional finite elements will be to 
investigate other methods of modelling three-dimensional bodies, to revise the 
Hrennikoff model to free the restriction on Poisson's ratio (v - l/4), and to 
develop Hrennikoff shapes other than the rectangular parallelepiped. 
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APPENDIX A 


3-D Hrennikoff Element Stiffness Matrix 


The stiffness matrix for the Hrennikoff element is derived in a straightfor- 
ward manner by considering the stiffness of each of the individual rods that make 
up the element and then assembling these into the overall stiffness matrix for the 
element. For the Hrennikoff element as shown in Figure 1 (Page 2) there are at 
most six different rods since all of those orientated parallel to any one axis are 
the same and all of those parallel to any plane are the same as indicated by the 
area equations in Section 2. 

Using the notation of Equations (2) the stiffness matrix for any rod in the 
Hrennikoff element written in a local coordinate system with axis along this rod 
can be written as: 


k 


p 




A-(l) 


where p takes on the values i , j , and k to give the stiffness matrices for those 
rods parallel to either the x i , x., or x k axes (see Figure 1) and 


k 


mn 


A E 

mn 


l 

-l 



A-(2) 


where mn take on the values i j , j k, ai,.d k i to give the stiffness matrices for 
the rods parallel to the x £ -Xj, x.-x k , and x (t -x i planes. The unassembled 
stiffness matrix for the Hrennikoff element is a matrix with the rod stiffness 
matrices on its diagonal and zero's elsewhere, denoting this as K we have: 


K 




A-(3) 
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where 


K. 


A-(4) 



and k A is givenby Equation A-(l). Similar definitions hold for and K^. . For 
K i . , K. k , and K ki these submatrices are defined similar to A-(4) but with 
A- (2) giving the form of the four diagonal matrices. 

Equation A-(3) contains the stiffness matrices of all 24 rod elements (four 
each of six types) in coordinate systems oriented along each rod. In order to 
get the stiffness matrix for the Hrennikoff element in the , x. , and x k coordi- 
nates we need to perform a coordinate transformation and apply joint (grid point) 
compatibility conditions, that is, we need to relate the individual rod degrees of 
freedom inherent in Equations A-(l) and A-(2) in local coordinates to the grid 
point degrees of freedom in the basis coordinates x. , x., and x k . This is ac- 
complished by a pre and post multiplication of K in Equation A-(3) by a trans- 
formation matrix to give: 


K = C T K C 


A- (5) 


where K is the Hrennikoff element stiffness matrix (24 x 24) in terms of dis- 
placements in the x. , x. , and x k coordinates of Figure 1 and C is a transfor- 
mation matrix whose form is: 


c. 
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c. . 

1 J 

c. . 

1 J 
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C iU 


Sci 
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Cui 


The submatrices C. ,_C. etc. are 8x8 matrices that relate the degrees of free- 
dom inherent in K £ , K. etc. to the grid point degrees of freedom of the Hrennikoff 
element. These submatrices in Equation A-(6) will not be presented here but are 
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easily found by geometrical considerations dependent on the ordering of the rod 
degrees of freedom and the Hrennikoff element corner point (grid point) degrees 


K 


K. . 

1 1 


symm. 


we have: 


K. . 

i j 

K. . 

i k 
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j j 

K M, 

i k 
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where K. . = [ C T K. C. + CT. K. . C + C. T K. . C . 1 

li ill lj ij ij ki ki ki 


A- (8) 


K. . = [ C T K. . C 

IJ I 1 J IJ 1 J 
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A-(10) 


K. . = [ C T K. C. + C T . K. . C. + C.T K . C.l 

J J [ J 3 3 1 J 1 3 13 3^ J k j kl 


A-(ll) 


K., = C. T . K.. C.J 

j k ^ j k j k j k J 


A-(12) 


^ = [<k T \ C* * c h K , k C jk t C.T K,, c kl ] 


A-(13) 


The 8x8 matrices give the stiffness coefficients of the Hrennikoff element. For 
example, the coefficients in the K. . matrix are the stiffness coefficients giving 
the forces in the x i direction at all eight grid points due to unit displacements 
of each of the eight grid points in the x t direction. Using A-(7) the stiffness 
equation for the Hrennikoff element is: 
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A-(14) 


Where, for example, r . , and R £ are the displacement and force vector compo- 
nents in the x. direction for the eight grid points of the Hrennikoff element. We 
will take the following definitions for the r. , r . , and r k displacement vectors: 
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Where, for example, r hi indicates the displacement at grid point h (see Fig- 
ure 1) in the x. direction. Using this definition and Equations A-(l) through 
A-(3) and A-(8) through A-(13) and also A-(15) the rod areas given by Equations 
(2) and (3) of Section 2, the stiffness matrices in Equation A-(14) are found to 
be as shown in Table 1 . 


The terms a i . , a ki are defined as: 
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also 


V = a { a. a k = Volume contained by the 
Hrennikoff element. 
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Table 1 


Hrennikoff Element Stiffness Matrix 
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Table 1 (Continued) 
Hrennikoff Element Stiffness Matrix 
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Table 1 (Continued) 
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